/* How to read indices:
- idx_* always refers to a neuron index.")
- the suffix tells you what the index is relative to,")
  e.g. idx_sim means "index relative to the entire simulation, irrespective of which realization"")
- similarly for other suffixes (div=division, realization, model, etc.)")
- prefixes are explained in comments.
*/

// >>> SimulationUnit._auxiliary_kernels(g) <<<

// The step kernel is called by the runtime on each timestep.
__kernel void step(const int t, const int realization_start,
                   __global float * const LIF_v,
                   __global float * const LIF_t_release) {
    // >>> SimulationUnit._enter_step_kernel(g) <<<

    // = Loop bound calculations
    // Index of this device thread.
    const int gid = get_global_id(0);
    // Total number of threads in division.
    const int gsize = get_global_size(0);

    // Index, relative to the simulation, of the first neuron being processed by the current division.
    const int first_idx_sim = realization_start * 4000;
    // and the last one
    const int last_idx_sim = min(first_idx_sim + 4000, 4000);

    // Loops over all neurons in this division (irrespective of which realization they are in.
    for (int idx_sim = first_idx_sim + thread_num; idx_sim < last_idx_sim; idx_sim += gsize) {
        // = Neuron index calculations
        // index of the realization this neuron is a part of
        // (remember how integer division works)
        const int realization_num = idx_sim / 4000;
        // index of the first neuron in this realization, relative to the sim
        const int realization_start_idx_sim = realization_num * 4000;
        // index of the first neuron in this realization, relative to the division
        int realization_start_idx_div = (realization_num - realization_start) * 4000;
        // index of this neuron relative to the realization
        int idx_realization = idx_sim - realization_start_idx_sim;
        // index of this neuron relative to the division
        int idx_division = idx_sim - first_idx_sim;

        // = Determine which model to invoke and run its generated code.
        // MODEL LIF
        // >>> ModelUnit._enter_model_code(mg) <<<
        // = Model index calculations
        // Index of neuron relative to the model.
        int idx_model = idx_realization - 0;

        // Index into state variable array for this model.
        // (state variables are contiguously arranged for all realizations of
        // this model in this division.)
        int idx_state = idx_model + (realization_num - realization_start)*4000;

        // >>> ModelUnit._pre_read_state_variables(mg) <<<
        // >>> ModelUnit._read_state_variables(mg) <<<
        const float v = LIF_v[idx_state];
        const float t_release = LIF_t_release[idx_state];
        // >>> ModelUnit._post_read_state_variables(mg) <<<

        // >>> ModelUnit._pre_read_incoming_spikes(mg) <<<
        // >>> ModelUnit._read_incoming_spikes(mg) <<<
        // >>> ModelUnit._post_read_incoming_spikes(mg) <<<

        // >>> ModelUnit._pre_calculate_inputs(mg) <<<
        // >>> ModelUnit._calculate_inputs(mg) <<<
        // >>> ModelUnit._post_calculate_inputs(mg) <<<

        // >>> ModelUnit._pre_process_incoming_spikes(mg) <<<
        // >>> ModelUnit._process_incoming_spikes(mg) <<<
        // >>> ModelUnit._post_process_incoming_spikes(mg) <<<

        // >>> ModelUnit._pre_variable_calculations(mg) <<<
        // >>> ModelUnit._variable_calculations(mg) <<<
        float v_new = t*0.1 >= t_release ? v + 0.1/20.0*(-v + currents) : v;
        // >>> ModelUnit._post_variable_calculations(mg) <<<

        // >>> ModelUnit._pre_independent_variable_updates(mg) <<<
        // >>> ModelUnit._independent_variable_updates(mg) <<<
        // >>> ModelUnit._post_independent_variable_updates(mg) <<<

        // >>> ModelUnit._pre_spike_processing(mg) <<<
        if (v_new >= 10.0) {
            // >>> ModelUnit._enter_spike_generated(mg) <<<

            // >>> ModelUnit._pre_spike_variable_updates(mg) <<<
            // >>> ModelUnit._spike_variable_updates(mg) <<<
            LIF_v[idx_state] = 0.0;
            LIF_t_release[idx_state] = t_release + 5.0;
            // >>> ModelUnit._post_spike_variable_updates(mg) <<<

            // >>> ModelUnit._pre_spike_propagation(mg) <<<
            // >>> ModelUnit._spike_propagation(mg) <<<
            // >>> ModelUnit._post_spike_propagation(mg) <<<

            // >>> ModelUnit._exit_spike_generated(mg) <<<
        } else {
            // >>> ModelUnit._enter_no_spike_generated(mg) <<<

            // >>> ModelUnit._pre_no_spike_variable_updates(mg) <<<
            // >>> ModelUnit._no_spike_variable_updates(mg) <<<
            LIF_v[idx_state] = v_new;
            // >>> ModelUnit._post_no_spike_variable_updates(mg) <<<

            // >>> ModelUnit._exit_no_spike_generated(mg) <<<
        }
        // >>> ModelUnit._post_spike_processing(mg) <<<
        // >>> SimulationUnit._exit_step_kernel(g) <<<
    }
}


